tsi_1984-2021
Shawn Schneidereit
2022-08-18
## test visualizazion of GreenLand masks
link <- "/data/Dagobah/greengrass/grassland_ger/01_CTM_GL_mask_17-19/GL_mask_2017-2019.tif"
mask_GL <- terra::rast(link)
plot(mask_GL)TSI I/C/O
Import TSI images (single location)
# Full mosaic
#link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt.ovr"
# base case
#link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.vrt.ovr"
# Single tile base case
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.tif"
tsi_base_raw <- rast(link)
# Single tile case day_16_sigma_8_16_32
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
tsi_day_16_sigma_8_16_32_raw <- rast(link)
# Single tile case day_16_sigma_8_8_8_16_16_16_64
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_8_8_16_16_16_64/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
tsi_day_16_sigma_8_8_8_16_16_16_64_raw <- rast(link)
# Single tile case day_16_sigma_8_8_8_16_16_16_32_64
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_8_8_16_16_16_32_64/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
tsi_day_16_sigma_8_8_8_16_16_16_32_64_raw <- rast(link)
# combine TSIs into a list to facilitate functional loop wise data operations
tsi_list <- list(tsi_base_raw, tsi_day_16_sigma_8_16_32_raw, tsi_day_16_sigma_8_8_8_16_16_16_64_raw, tsi_day_16_sigma_8_8_8_16_16_16_32_64_raw)
cases <- c("base_case", "day_16_sigma_8_16_32", "day_16_sigma_8_8_8_16_16_16_64", "day_16_sigma_8_8_8_16_16_16_32_64")
names(tsi_list) <- cases
tsi_obs_length <- list()
lapply(seq_along(tsi_list),
function(i)
tsi_obs_length[i] <<- length(names(tsi_list[[i]]))
)test vis of TSI
create stack for direct location comparison
tsi_stack <- c(tsi_list[[1]][[c(430)]],
tsi_list[[2]][[c(675)]],
tsi_list[[3]][[c(675)]],
tsi_list[[4]][[c(675)]])static maps
#de <- st_read("/data/Dagobah/dc/deu/vector/DEU.gpkg")
#plot(tsi_list[[2]][[800:808]])
lapply(seq_along(tsi_list),
function(i) plot(tsi_stack[[i]],
legend=T, col = viridis(n=100,option="D"),
range=c(0,10000),
plg=list(title= paste(names(tsi_list[i]), "\n", names(tsi_stack[[i]]))
))
)
Interactive map
pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(tsi_list[[4]][[675]]),
na.color = "transparent")
leaflet() |>
addTiles() |>
addRasterImage(tsi_list[[4]][[675]], group = "show_layer") |>
addLegend(pal = pal,
values = values(tsi_list[[4]][[600:600]]),
title = "Reflectance") |>
addLayersControl(
overlayGroups = c("show_layer"),
options = layersControlOptions(collapsed = FALSE)
)data wrangeling of TSI images
Compute and set min and max values for tsi cases
lapply(seq_along(tsi_list),
function(i)
tsi_list[[i]] <<- setMinMax(tsi_list[[i]])
)## [[1]]
## class : SpatRaster
## dimensions : 1000, 1000, 675 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 4166026, 4196026, 3344920, 3374920 (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
## source : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.tif
## names : 19840~LND05, 19840~LND05, 19840~LND05, 19840~LND05, 19841~LND05, 19841~LND05, ...
## min values : -507, 1140, 2569, 1768, 716, -355, ...
## max values : 8815, 8873, 8598, 8288, 7954, 9881, ...
##
## [[2]]
## class : SpatRaster
## dimensions : 1000, 1000, 866 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 4166026, 4196026, 3344920, 3374920 (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
## source : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif
## names : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ...
## min values : NaN, NaN, NaN, -507, 475, 225, ...
## max values : NaN, NaN, NaN, 8815, 8873, 8873, ...
##
## [[3]]
## class : SpatRaster
## dimensions : 1000, 1000, 866 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 4166026, 4196026, 3344920, 3374920 (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
## source : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif
## names : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ...
## min values : 570, 592, 613, 871, 875, 349, ...
## max values : 8873, 8873, 8873, 8873, 8873, 8873, ...
##
## [[4]]
## class : SpatRaster
## dimensions : 1000, 1000, 866 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 4166026, 4196026, 3344920, 3374920 (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
## source : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif
## names : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ...
## min values : 570, 592, 613, 703, 857, 448, ...
## max values : 8873, 8873, 8873, 8873, 8873, 8873, ...### alternatice slower na layer computation method
#
# na_index <- list("base_case"=c(), "day_16_sigma_8_16_32"=c(), "day_16_sigma_4_8_16"=c(), "day_16_sigma_8_8_8_16_16_16_64"=c(), "day_8_sigma_8_16_32"=c())
#
# system.time(
#
# lapply(seq_along(tsi_list),
# function(i)
#
# for (j in c(1:nlyr(tsi_list[[i]]))) {
#
# na_index[[i]][[j]] <- minmax(is.na(foo[[i]][[j]]))[1]
#
# }
# )
#
# )For all cases most cases, there are occasionally layers consisting exclusively of NAs (ie where min and max is NA). Overall there are relatively rare (~3-15%; computed below), and occur mostly at the beginning of the Landsat program (1984-1986), and occasionally during winter months.
The empty layers are included in the force output to satisfy the specified day window intervals for the interpolation. Unfortunately this has the downstream effect that when attempting to convert the spatial object to a df, R throws errors when empty layers are included. The current workaround is to compute the the min/max per layer, create an index of valid layers, and then use this index to subset the tsi stacks.
remove all NA layers
# create filter out all layers that dont contain any values
tsi_na_index <- list()
lapply(seq_along(tsi_list),
function(i)
tsi_na_index[i] <<- minmax(tsi_list[[i]]) |>
as.data.frame() |>
slice(1) |>
pivot_longer(everything()) |>
mutate(id = row_number()) |>
na.omit()
)
# convert into index vector
lapply(seq_along(tsi_list),
function(i)
tsi_na_index[[i]] <<- as.vector(tsi_na_index[[i]])
)
# subset tsi object with only layers containing values
lapply(seq_along(tsi_list),
function(i)
tsi_list[[i]] <<- tsi_list[[i]][[ tsi_na_index[[i]] ]]
)
# filter out reflectance values below and above 0-100 % range
# lapply(seq_along(tsi_list),
# function(i)
# tsi_list[[i]] <- tsi_list[[i]] |>
# filter(across(1:length(names(tsi_list[[i]])), ~ . > 0),
# across(1:length(names(tsi_list[[i]])), ~ . < 10000))
# )Check NA layer occurance
# lapply(seq_along(tsi_list),
# function(i)
# tsi_obs_length[[i]] - length(names(tsi_list[[i]]))
# )
lapply(seq_along(tsi_list),
function(i)
(paste("In the", names(tsi_list)[[i]],"scenario", (1-round(( length(names(tsi_list[[i]]))/tsi_obs_length[[i]]), digits=4))*100, "% of layers consisted completely of NA values (n=", ( tsi_obs_length[[i]] - length(names(tsi_list[[i]]))), "out of", tsi_obs_length[[i]], ")"))
)[[1]] [1] “In the base_case scenario 0.149999999999995 % of layers consisted completely of NA values (n= 1 out of 675 )”
[[2]] [1] “In the day_16_sigma_8_16_32 scenario 3.35 % of layers consisted completely of NA values (n= 29 out of 866 )”
[[3]] [1] “In the day_16_sigma_8_8_8_16_16_16_64 scenario 0 % of layers consisted completely of NA values (n= 0 out of 866 )”
[[4]] [1] “In the day_16_sigma_8_8_8_16_16_16_32_64 scenario 0 % of layers consisted completely of NA values (n= 0 out of 866 )”
create updated stack for direct location comparison
(ie because NA filter removed some layers and dates dont correspond when using the previous index)
tsi_stack <- c(tsi_list[[1]][[c(430)]],
tsi_list[[2]][[c(648)]],
tsi_list[[3]][[c(675)]],
tsi_list[[4]][[c(675)]])#rasterVis::levelplot(tsi[[c(24)]], par.settings = viridisTheme)
#rasterVis::levelplot(tsi[[c(1,4,24,25:35)]], par.settings = viridisTheme)
lapply(seq_along(tsi_list),
function(i)
plot(tsi_stack[[i]],
legend=T, col = viridis(n=100,option="D"),
range=c(0,10000),
plg=list(title= paste(names(tsi_list[i]), "\n", names(tsi_list[[i]][[674]]) ))
)
)generation of sample points (performed in Rstudio GUI)
roi_names <- c("watt_GL", "large", "medium", "small")
# select the ROI for each location type
# watt_GL <- terra::sel(tsi_list[[4]][[c(430:430)]])
# large <- terra::sel(tsi_list[[4]][[c(430:430)]])
# medium <- terra::sel(tsi_list[[4]][[c(430:430)]])
# small <- terra::sel(tsi_list[[4]][[c(430:430)]])
# used to make raster stack (c() stack does not work due to extents not matching)
#roi_list <- c(watt_GL, large, medium,small)
#save roi sample points raster stack
# writeRaster(watt_GL, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/watt_GL.tif",
# filetype='GTiff', overwrite=TRUE)
# writeRaster(large, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/large.tif",
# filetype='GTiff', overwrite=TRUE)
# writeRaster(medium, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/medium.tif",
# filetype='GTiff', overwrite=TRUE)
# writeRaster(small, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/small.tif",
# filetype='GTiff', overwrite=TRUE)
#load sampled points
watt_GL <- rast("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/watt_GL.tif")
large <- rast("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/large.tif")
medium <- rast("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/medium.tif")
small <- rast("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/small.tif")
# used to make raster stack (c() stack does not work due to extents not matching)
#roi_list <- c(watt_GL, large, medium,small)
roi_list <- list(watt_GL, large, medium,small)
# Convert each ROI to df
#
# watt_GL_df <- watt_GL |> as.data.frame(xy=T) |> select(x,y) |> mutate(type="watt_GL")
# large_df <- large |> as.data.frame(xy=T) |> select(x,y) |> mutate(type="large")
# medium_df <- medium |> as.data.frame(xy=T) |> select(x,y) |> mutate(type="medium")
# small_df <- small |> as.data.frame(xy=T) |> select(x,y) |> mutate(type="small")
# Combine locations into single df
# sampled_points <- rbind(watt_GL_df, large_df, medium_df, small_df)
#save sampled points
#saveRDS(sampled_points, file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/sampled_points.rds")
# load sampled points
sampled_points <- readRDS("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/sampled_points.rds")Plots of ROIs
sampled location (little red highlighed areas)
knitr::include_graphics("~/Documents/R/HU_SHK/tsi_1984-2021/ROI_location_map.png")
Detailed plots of each ROI
lapply(seq_along(roi_list),
function(i)
plot(tsi_list[[4]][[c(430:430)]],
legend=T,
col = viridis(n=100,option="D"),
range=c(0,10000),
ext=ext(roi_list[[i]])*10,
plg=list(title= paste(roi_names[i], "\n", names(tsi_list[[i]][[600:600]]))
))
)
compairison of TSI cases for each location
lapply(seq_along(roi_list),
function(i)
plot(tsi_stack,
legend=T,
col = viridis(n=100,option="D"),
ext=(ext(roi_list[[i]])*10),
plg=list(title= paste(names(tsi_list[i]), "\n", roi_names[i])),
range=c(0,10000))
)
For the displayed date (July 2013), all plots seem to show quite congruous results. As expected the un-interpolated time series shows more data gaps, but in general the visual correspondence of reflectance values look good.
TSI visualization
data manipulation and subsample creation
Due to C memory usage limit issues, the conversion from spatraster to spatial data frame had to be split into small individual processing chunks, that are then subsequently joined to form the single df. First unique start and stop limits for the processing loop are defined (as cases have heterogenious lengths depending on number of NA layers and date window size). Following this the range of the processing checks is set in a list. Then the spat raster object of each case are converted into spatial df´s, via the defined subset processing chunks. Finally, all the individual chunks of a case are joined into a single spatial df object
… works but is not a pretty workaround…
tsi_df <- list("base_case"=c(), "day_16_sigma_8_16_32"=c(), "day_16_sigma_4_8_16"=c(), "day_16_sigma_8_8_8_16_16_16_64"=c(), "day_8_sigma_8_16_32"=c())
start <- list()
stop <- list()
for (i in 1:length(names(tsi_list))) {
# create processing chunks
start[[i]] <- seq(from = 1, to = length(names(tsi_list[[i]]))-20, by = 20)
stop[[i]] <- seq(from = 20, to = length(names(tsi_list[[i]])), by = 20)
}
# create a unique length of processing chuncks based on the number of observations in each scenario
processing_chunks <- start[[1]][c(1:length(start))] # initial list
# loop though each case
for (j in 1:length(start)) {
# loop through each start and stop combination
for (i in 1:length(start[[j]])) {
processing_chunks[[j]][i] <- list(c(start[[j]][i], stop[[j]][i]) )
}
}
# Convert spat raster into df, via subset processing chunks
### uncomment to execute
# lapply(seq_along(tsi_list),
# function(i)
#
# for (j in c(1:length(processing_chunks[[i]]))) {
#
# tsi_df[[i]][[j]] <<- terra::as.data.frame(tsi_list[[i]][[c(processing_chunks[[i]][[j]][1]:
# processing_chunks[[i]][[j]][2])]],
# xy = TRUE, na.rm=F) |>
# # filter out rows that only contain na
# filter_at(vars(-x,-y), any_vars(! is.na(.))) |>
#
# left_join(sampled_points) |>
# filter(!is.na(type))
# }
# )
# create join for list into df
# tsi_df_joined <- list()
# lapply(seq_along(tsi_df),
# function(i)
# tsi_df_joined[[i]] <<- reduce(tsi_df[[i]], full_join) |>
# group_by(type) |>
# sample_n(1)
# )
#save tsi df
#saveRDS(tsi_df_joined, file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/tsi_df_joined.rds")
# load tsi df
tsi_df_joined <- readRDS("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_df_joined.rds")
# create long version of df with unique and formatted location and date cols
tsi_df_long <- list()
lapply(seq_along(tsi_list),
function(i)
tsi_df_long[[i]] <<- tsi_df_joined[[i]] |>
pivot_longer(cols = c(-x,-y,-type), names_to = "timestamp") |>
mutate(xy= paste0(x,"_",y)) |>
mutate(date = as.Date(timestamp, format="%Y%m%d"),
year = lubridate::year(date),
value = if_else(value<0,0,as.numeric(value)),
case = names(tsi_list[i]))
)plot faceted time series 1984:2000
lapply(seq_along(tsi_list),
function(i)
tsi_df_long[[i]] %>%
filter(year<2001) %>%
ggplot( aes(date, value, color=type, group=xy)) +
geom_line(aes(group=type)) +
geom_point(size=0.5) +
theme_classic() +
scale_colour_viridis_d(option = "D") +
scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust=1),
legend.position="bottom"
) +
ggtitle(cases[i]) +
facet_wrap(~year, scales = "free_x")
)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
plot faceted time series 2001:2021
lapply(seq_along(tsi_list),
function(i)
tsi_df_long[[i]] %>%
filter(year>2001) %>%
ggplot( aes(date, value, color=type, group=xy)) +
geom_line(aes(group=type)) +
geom_point(size=0.5) +
theme_classic() +
scale_colour_viridis_d(option = "D") +
theme_minimal() +
scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme(axis.text.x = element_text(angle = 45, hjust=1),
legend.position="bottom"
) +
ggtitle(cases[i]) +
facet_wrap(~year, scales = "free_x")
)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
lapply(seq_along(tsi_list),
function(i)
ggplot(tsi_df_long[[i]], aes(date, value)) +
geom_smooth() +
theme_minimal() +
scale_x_date(date_breaks = "1 year", date_labels = "%d-%m-%y")+
theme(axis.text.x = element_text(angle = 45 , hjust=1))
)lapply(seq_along(tsi_list),
function(i)
tsi_df_long[[i]] %>%
ggplot( aes(date, value, group = year)) +
geom_boxplot() +
theme_minimal() +
scale_x_date( date_labels = "%d-%m-%y")+
theme(axis.text.x = element_text(angle = 45, hjust=1))
)direct location and case compairison
tsi_df_long_combined <- reduce(tsi_df_long, full_join)## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")
## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")
## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")ggplot(tsi_df_long_combined, aes(date, value, color=case, group = xy)) +
geom_line(aes(group=case), alpha=0.5) +
geom_point(size=0.5, alpha=0.9) +
geom_smooth(color="black") +
theme_minimal() +
scale_colour_viridis_d(option = "D") +
scale_x_date(date_breaks = "4 year", date_labels = "%Y")+
theme(axis.text.x = element_text(angle = 45, hjust=1),
legend.position="bottom") +
# ggtitle(cases[i]) +
facet_grid(rows = vars(type), cols = vars(case), scales="free_y") ## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#facet_grid( ~year)
#facet_grid(vars(drv), vars(cyl))compairson of Base case values and interpolated TSI outputs
year_chunks <- round(seq(1984, 2021, length.out = 5))
lapply(seq_along(c(1:4)),
function(i)
tsi_df_long_combined %>%
filter(year > year_chunks[i],
year < year_chunks[i+1],
case != "base_case") %>%
ggplot(aes(date, value, color=case, group = xy)) +
geom_line(aes(group=case)) +
geom_point(size=0.5) +
geom_point(data=tsi_df_long_combined |> filter(case =="base_case",
year > year_chunks[i],
year < year_chunks[i+1]),
aes(date, value),
color ="red", alpha=0.5) +
theme_minimal() +
scale_colour_viridis_d(option = "D") +
scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme(axis.text.x = element_text(angle = 45, hjust=1),
legend.position="bottom") +
facet_grid(rows = vars(type), cols = vars(year), scales="free_x")
# facet_wrap(~year, scales = "free_x")
)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
plot time series with formatted dates
ggplot(tsi_df_long[[i]], aes(date, value, group=xy, color=xy)) +
geom_line() +
theme_classic() +
scale_colour_viridis_d(option = "D") +
scale_x_date(date_breaks = "1 month", date_labels = "%m-%y")+
theme(axis.text.x = element_text(angle = 45, hjust=1),
legend.position="none")
ggplot(tsi_df_long[[i]], aes(date, value)) +
geom_smooth() +
theme_minimal() +
scale_x_date(date_breaks = "1 year", date_labels = "%d-%m-%y")+
theme(axis.text.x = element_text(angle = 45 , hjust=1))
ggplot(tsi_df_long[[i]], aes(date, value, group = timestamp)) +
geom_boxplot() +
theme_minimal() +
scale_x_date( date_labels = "%d-%m-%y")+
theme(axis.text.x = element_text(angle = 45, hjust=1)) ggplot + tidyterra vis option (nice but runs slow)
ggplot() +
geom_spatraster(data = tsi[[c(4:8)]]) +
facet_wrap(~lyr, ncol = 2) +
scale_fill_whitebox_c(
palette = "viridi",
# labels = scales::label_number(suffix = "º")
) +
labs(fill = "NDVI") +
theme_classic()## gif image creation
for (i in c(1:100
#length(names(tsi))
)) {
# Step 1: Call the pdf command to start the plot
png(file = paste0("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/ts_gif/day_16_sigma_8_16_32/plot_", names(tsi_list[[1]])[[i]], ".png"), width = 265, height = 125, units='mm', res = 100) #, # The directory you want to save the file in
# width = 12, # The width of the plot in inches
#height = 12) # The height of the plot in inches
plot(tsi_list[[2]], legend=T, col = viridis(n=100,option="D"),
range=c(0,10000),
plg=list( title = paste(as.Date(names(tsi_list[[1]][[i]]), format="%Y%m%d")), title.cex=1.25)
)
# Step 3: Run dev.off() to create the file
dev.off()
}TSI gif
knitr::include_graphics("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/ts_gif/day_16_sigma_8_16_32/ts_16_8_16_32.gif")
# old attempted at looped df processing
processing_chunks <- list(c(1,10),c(11,20),c(21,30),c(31,40),c(41,50))
processing_chunks <- c(1,100,200,300,400,500,600,700,800)
processing_chunks <- c(1,10,20,30) # short version for testing
tsi_df <- list("base_case"=c(), "day_16_sigma_8_16_32"=c(), "day_16_sigma_4_8_16"=c(), "day_8_sigma_8_16_32"=c())
for (j in c(1:3
#(length(processing_chunks)-1)
)) {
lapply(seq_along(tsi_list),
function(i)
tsi_df[[i]][[j]] <<- terra::as.data.frame(tsi_list[[i]][[c(processing_chunks[[j]][1]:
processing_chunks[[j]][2])]],
xy = TRUE, na.rm=F) |>
# filter out rows that only contain na
filter_at(vars(-x,-y), any_vars(! is.na(.))) |>
# random subset for data vis
slice_sample(n=15)
)
}
foo <- list()
lapply(seq_along(tsi_df),
function(i)
foo[[i]] <<- reduce(tsi_df[[i]], left_join)
)### test for R tile merging
link_a <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
link_b <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0054_Y0052/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
a <- rast(link_a)
b <- rast(link_b)
ab <- merge(a,b)
plot(ab[[c(10:13)]], legend=F, col = viridis(n=100,option="D"),
range=c(0,10000)
)
object.size((ab))
writeRaster(ab, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/test_mosaic_r.tif",
# overwrite=TRUE,
datatype = "INT4S",
filetype='GTiff')
link_ab <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/test_mosaic_r.tif"
ab <- rast(link_ab)
plot(ab[[c(10:13)]], legend=F, col = viridis(n=100,option="D"),
range=c(0,10000)
)